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ABSTRACT 

The purpose of this paper is to establish a method for identifying un- 
known parameters involved in the boundary state of a class of diffusion sys- 
tems under noisy observations. A mathematical model of the system dynamics is 
given by a two-dimensional diffusion equation, whose boundary condition is 
partly unknown due to the existence of an unknown parameter. Noisy observa- 
tions are made by sensors allocated on the system boundary. Starting with the 
mathematical model mentioned above, an on-line parameter estimation algorithm 
is proposed within the framework of the maximum likelihood estimation. Exis- 
tence of the optimal solution and related necessary conditions are dis- 
cussed. By solving a local variation of the cost functional with respect to 
the perturbation of parameters, the estimation mechanism is proposed in a form 
of recursive computations. Finally, the feasibility of the estimator proposed 
here is demonstrated through results of digital simulation experiments. 


Research was supported under the National Aeronautics and Space Administration 
under NASA Contract No. NAS1-18107 while the second author was in residence at 
ICASE, NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 


Recently, there has been much Interest In the parameter identification 
problem for distributed parameter systems. For parameter identification prob- 
lems for a class of distributed parameter systems, excellent survey papers 
have been published by many researchers (e.g., Kubrusly, 1976). Among them, 
the method for parameter estimation including discussions on convergence 
properties of computer algorithm has been developed by Banks and Lamm (1985), 
Kravaris and Seinfeld (1985), etc. Recently, feasible methods for estimating 
unknown coefficients which depend on the state and spatial variables were re- 
ported by the authors. In this paper, motivated by parameter identification 
problems of such parameters as heat conduction in chemical reactor or oil 
reservoir problems (e.*., Chavent, 1978), we deal with the identification of 
boundary parameters for a two-dimensional diffusion system under noisy obser- 
vations . 


2. PROBLEM CONSIDERED 

bet T and G be the time interval [0,tf] and system domain given by a 
o 

bounded set in R , and let 8C be the boundary of domain C. The bound- 
ary 30 is characterized by C^-parameterized curve decomposed into two 
parts such that 

3G - 3G U 3G . (2.1) 

U K 

The system state u(t,x) at time t and at the location x = (xj,X2) is 
governed by heat diffusion equation: 


* 
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3U 3t >X> ~ Au ^» x > " 0 in T xG (2.2a) 

together with the Initial and boundary conditions 

u(0,x) - Uq(x) on G, (2.2b) 

+ O - 9(0)u(t,0 - o (2.2c) 

on T x aGy 

~ a ' n ~ = g(C) on T x 3 G r (2. 2d) 

2 7 2 2 

where A = 9 /3x‘ + 9 9 /9 n denotes the outer normal derivative and 

0(C) is unknown parameter with a form of 

0(0 = 1 /(b(C) + 1) (2.3) 

and where b(C ) denotes the heat transfer coefficient. Since the heat 
transfer coefficient b has its physical value of positive definite, then 
0 < 0(C) < 1. It is apparent that, in the case where 0(C) ■ 0, the bound- 
ary condition becomes the Dlrichlet type while, in the case where 0(C) * 
we have the boundary condition of Neumann type. 

Observation data are acquired through the sensors located on the boundary 

part 9G , as shown in Fig. 1. The observation mechanism is modeled by 
K. 
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O ^Allocation of Sensors 


Fig. 1. Location of sensors on the boundary 


t 

y(t) = / H[u(s; 0 ) ]ds + v(t), (2.4) 

0 

where y(t) - fy,(t),*»«,y (t)]' and v(t) denotes the observation noise 

I IT! 

term modeled by an m-d imensional brownian motion process. In Eq . (2.4), H{ • ] 
denotes the signal process defined by 


H[u(t ,0)1 = [h lU (t,P K (l) ;9),“' 


h u(t, P <m) ;0)]' 
in iv 


(2.5) 


where h^ (I = m) are assumed to be some constant, and where 

denotes the coordinates of sensor locations (zj*\z^^). 


P 


(i) 

K 
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The 

problem cons ide red 

here Is to 

f Lnd a method 

for identifying 

parameter 

9(0 

defined 

on 3 Gy 

from information 

of the a priori 

boundary 

known input 

gOO 

and the noisy 

observation data 

{ y(t)} t >o- 


3. THEORETICAL RESULTS FOR PARAMETER ESTIMATION 

Let Uj and ^ be measures Induced on the processes y(t) and 
y(t) = v(t), respectively. Then, the Radon-Nikodym derivative Is 

dp , t 

- — = exp [ f H[u(s;6 ) ] 'dy(s) 

M 2 0 ° 

(3.1) 

, t 

—*■ J H[u(s;6 )]'H[u(s;0 )]ds] 
i J n o o 


where 0 ^ denotes the true function of unknown parameter. Associated with 

(3.1), the likelihood ratio functional is given by 


t 

A (u,v,9 ) = exp [/ H[u(s ;0 ) ]dy(s) 

0 


-y / H[u(s; 0 ) ] "H [u(s ;0 ) Ids] . 
0 


(3.2) 


The maximum likelihood estimate for the parameter 0 is the solution of 


Inf 1(6) = I (0°) 
0 e H C 


where (H) denotes the admissible parameter class and 


(3.3) 
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1^(9 ) “ ~ — lnA t (u,y,0). 


(3.4) 


(See Baladrishnan, 1975, for more details.) By using the gradient method (see 
Polak., 1971), the optimal solution can be obtained by the following recursive 
computations: 

9 (l+1) = + X.V 0 I t (9 (i) ) (3.5a) 


e (0) - e 


for i=l , 2, • • • , 


(3.5b) 


where V denotes the gradient of the cost I (0). In order to compute 
0 t 

Vgl (0^^), we have to solve the following partial differential equation for 


each 0 


with 


(t). 


- Au(t ,x) =0 in T x G 


u(0,x) = Uq(x; 0^^) on G 


(3.6a) 


(3.6b) 


e (1) (f) + O _ f * (i) ( p ,))u(t,0 - o 


on T x 3G, 


9u(t , g) 
3n 


U 


g(5) on T x 9G r . 


(3.6c) 


(3.6d) 


This implies than in order to solve (3.3) massive computations are required. 
A possible method for avoiding this difficulty is to replace H[u(t;0)] in 
(3.2) bv the stationary value H [TT( 0)1 where TT(x;0) * lim u(t,x;0) and 
where u(x;0) is a solution of 


Au(x) * 0 in C 


(3.7a) 
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with boundary conditions 

0 (5 )~^- + (1 - e(5))u(C) - o on 3G U 

^i=g(£) on 3G k . 

Thus, we introduce the following cost functional: 

J t (0) = - 1 H[u(0)]'y(t) 

+ j H [ u( 0 ) ] ' H{u(8)]. 

The following result states the relation between the likelihood ratio 
t tonal and the cost functional, in this paper. 

Lemma Is (See Sunahara, et al. , 1984) The cost functional (3.8) 
fies the following asymptotic property : 

Lim E{ |J (0) - I t (0)| 2 } = 0. 

£**•«> 

* 

Thus, our problem is to find the optimal solution 0 satisfying 

inf J (0) = J (0*). 

oe(|t> 


(3.7b) 

(3.7c) 

(3.8) 
func- 

satls- 

(3.9) 

(3.10) 


In this paper, we shall define the admissible class of parameters as follows: 
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[Definition] Let (H) be the set of all admissible parameters. Then, (ff) is 
said to be the admissible parameter, i.e., 0 € (J), if 0 satisfies the 
following properties: 


(P-1) 


0 < 0 ( 5 ) < 0 < 1 on 9G, 


(P-2) 


8 € (T( aG 0 ). 


The next results give the existence property for the optimal solution of 
(3.10). 


Theorem 1: Suppose that 


(C-l) 


g € L (3G R ). 


Then, there exists at least one solution of (3.10) for a fixed t T with 


probabi 1 i ty one. 


The proof of Theorem 1 will be shown in Appendix 1. 


The next results show the necessary condition for the optimality of this 
parameter estimation problem. 


Theorem 2: Let 9 be the optimal solution of (3.10). Then the 

necessary condition for optimality is characterized by the following varia- 


1 


/ (0(0 - 0*(O) — 

3G (J {9 t (0-l} 


3u(C;9*) 3p(?;9*) 

—in — — rs — d? 

for V 0 , 9* C 9 


(3.11) 


inhere uO^) 

whe re pO^) 

nonhomogeneous 


is the solution of (3,7) corresponding to 9*0* and 
Is the solution of the following Laplace equation with the 
boundary condition, 


Ap(x) = 0 in G 


(3.12a) 


9*(5)i2|L ) -+ (1 - 9*(0)p(0 = 0 on 3 Gy 

*£1S2 = -n*[yi?l - H[u(8*)]) on 3G k . 


The proof of Theorem 2 will he shown in Appendix 2. 


(3.12b) 


(3.12c) 


By using the projection gradient method with the fixed step size A, the 
computational algorithm for solving (3.10) can be considered as follows: 


0 (l+1) = H 0 (0 (i) - A7 0 J t (0 (i) )) (3.13a) 

0 (0) = F for i = 1,2, 


• • • 


(3.13b) 
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where 11^ denotes the prolection operator on (U>, Then, applying Theorem 1, 


(3,13a) becomes 




TO 


3n 


9n 


(3.13c) 


In the sequel, we propose the on-line parameter estimator. Let t^ and t^ 
be the initial and the terminal time for parameter estimation. We choose the 
time step k as 


fc f “ C 0 


(3.14) 


and the time Interval [tf^tf] is divided into 


{t< n) >: 0 < e 0 - t<"> < t' n) < 


'• < tj n> < 


, (n) 

'• < t - t, 
n f 


* t 0 + ik for i » l,2,»»«,n- 


1 . 


(3.15) 


(n) 


by using the 


He compute the recursive algorithm (3.13) at each time t 
observed data y(t^ n b. Considering the fixed step size \ in (3.13) as 


the time step k, the estimated sequence 6 

.(n) 


:(t) 


can be obtained at each 


t ime t , 


for i = 1,2,»*» ,n. 


Theorem 3: The estimated sequence {8^0}? is characterized by the 

following variational Inequal Ity: 


f I 


for i * 1, 2,» • • ,n 


(3.16a) 
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where 


/ (0 (1+l) (C) - 0 (i) (O) 

3G,, 


x (<KO - 0 (1+l) (C))d? 


+ k r . , 1 ».“<€;P (i) > 

3G fJ {0 (i) (5)-l} 2 3n 
* (♦(?) - 0 (1+1) (O)dC >0 


for v e <M> 

i = 1, 2,» • • ,n-l 


(3.16b) 


i and p are the solutions of 


~(i) 

/ u(x)(Ai|/(x))dx = / — - 7 - j - v — 

G 3 Gy {0^ )~1> 


„ 3u(0 3<KO 

x 5n §n d5 


r 3uTc ) , Xjl . 
-/ TTT"* ( °d 5 


3G 


U 


*S u(5)^dc 

3C K 


- / g(5 )♦({)« 

d G 

K 

^ (i) 

{ p(x)(Ai|j(x))dx = / ~ ^Y — 

G 3Gy {0 U; (O“1} 


(3.17) 


3p(Q 3 »(S) 
8n 9n 




x 
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-j Mil ,(,« 

dG u 

3G k 

+ Hi*]' ( -- y(t t (n) ) - Hru(0 (i) )]) 
r i 

for V, e H 2 (G). 

♦ 


A. ESTIMATION ALGORITHM 
4.1. B o a ndary Element Approach 

In order to implement the proposed algorithm in the previous section, we 
adopt the boundary element me thod (BEM) (see Brebbia, 1978). Let v(x,x*) be 
the fundamental solution satisfying 

Avfx.x*) + 6(x - x*) ■ 0, (4.1) 

where 6(x - x*) is the Dirac delta function. It is well-known that 

v(x,x*) is given by 

v(x,x*) * In , (4.2) 

r(x,x ) 

where r(x,x*) denotes the distance between x * (x^,X 2 ) a °d x* * (xj^). 
By using Green's formula and from (4.2), Eq. (3.7) can be rewritten by the 
following boundary integral equations: 
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2it C u(x' ;0 ) 


+ / { 9TO=T IJ <ln 5 > 

3C„ 3n rtx 1 ,?) 


- ln _ ■ _) ’”«•»> d ? 

rC* 1 ,?) 3n 


+ / x— (In ■ ■ ? — )u(g ;9)dS 
3C„ 3n rCx 1 ,?) 


- / ln j— g ({)d5, 


3G k r(x ,£) 


where denotes the constant such that 


1 for x € G 


C. = < 1/2 for x € 9G 


0 for x € (T- . 


Similarly, from Eqs. (3.12), we obtain 


2ir C j ,p(x ;0 ) 


+ f f 9 ^> 9 _ fin 1 N 

+ J 77 'i- 1 in -r ' 


-1 9n , i _ • 

r(x ,5) 


- In { } <JC 

rC* 1 ,?) 3n 


+ / J- On J )p(C;9)d5 


r(x ,5) 
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-/ In 1 H*[^- - H[u(9)]]dC. (4.5) 

3G k rCx 1 ,£) 

In order to perform BEM, the boundary 3G is approximated by 3G, which 
is decomposed into (n + m) elements, i.e., as illustrated in Fig. 2 

SG-S'gylja^ (4.6a) 

3Gy - SG^U .*.U (4.6b) 

3G r - 3Gg 1} U ••• U 3G^ m ^ (4.6c) 



Fig. 2 Boundary elements and nodes of the system. 


In Fig. 2, { ^ j » i and denote the node of boundary elements 

{3Gy^}j m j and {3 g£^}” m j respectively. It is noted that the nodes 
{p£^}j_l coincide with the sensor allocation as shown in Fig. 1. 


For 
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For economy of notations, Introduce the following vector forms: 

{<KV } = I<l>(p5 1) ),***,<KPj n) )]' (4.7a) 

{4> (P K - N>(p£ 1) ),***,<Kp£ m) )K. (4.7b) 


Then, from (4.4) and (4.5), the boundary steady state Is comovit<-<1 by the fol- 
lowing vector-matrix equations: 


_{ u(P K ;e)} 


= A -1 (9 )C 


{ 0 } 

{i( p K ) } 


( l£ <V 9 » 

{ p(P K ;0)} 


{ 0 } 

- H[u(9)]]} 


where A(0) 


and C are (n + m) x (n + m) matrices. 


(4.8) 


(4.9) 


4.2 Estimation Algorithm 

The numerical procedure for the estimation algorithm in the previous sec- 
tion can be presented as follows: 


Step 0: Select an Initial parameter 



at time t - t^ > 0 such that 


S<°> 


{0 (O) (P U )} € Cfl), 


and set 


i 


0 
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.:d) 


Step 1: Compute {u(P u ;0 vw )} and (p^jO^)} by solving (4.8) and 

(4.9). 


* ( i +1 ^ 

Step 2: Compute the estimated parameter 0 V ' at time t * t Q + (i + l)k 


by 




+ X 


1 


{0(P^ J) ) - l 2 } 


au(P^ j) ;e) 3p(P l ( J j) ;0) 


9 n 


8n 


0=0^^ 


and 


for P*' 1 ' on 9 Gy 


; (i+ 1 ) <p<i>> - v < 1 +n)(p u 3)) - 


(4.11a) 


(4.11b) 


Step 3: Setting i by i + 1, go back to Step 1. 


5. NUMERICAL EXPERIMENTS 

Throughout numerical experiments, the system domain is given by a rec- 
tangle as illustrated in Fig. 3. Figure 3 shows also that the boundary is 
divided into 20 elements, while 9Gy and 3G R are decomposed into 5 and 

15 elements, i.e., n = 5 and m = 15, respectively. The boundary input g 
is set as 

on 9G^. 


g(0 - 500 


(5.1) 
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and where 


* (i) 
Av . 

1 


process generated 


state u( t ,x;a ) 


denotes an Increment of the standard Brownian 
by the Gaussian distribution. In Eq . (5.2) , the 

is computed by adopting BEM. 


motion 

system 


Example i: The unknown boundary parameter is preassigned as 0^ * 0.2* By 

substituting known boundary input g and observation data 

(i = l,2,«**) into the proposed algorithm, starting with an initial 

parameter at time tQ = 200 x k (k * 1.0), the estimated parameter 
{0 (E)}._a i ## could be computed. Results of this example are demon- 

strated in Fig. 4. 

e t (e>. a 0 (e> 



e 

Fig. 4. The >'~iLimar.ed parameter in Example 1. 

Example 2: In Example 2, we deal with the case where the unknown parameter is 

a space-variable function. The unknown parameter and simulation results are 
illustrated in Fig. 5. 
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fl t (e). e c (f ) 



rljr* *>. The estimated parameter in Example 2. 


6. CONCLUSIONS 

A method for identifying the boundary condition has been developed in 
this paper for a distributed system whose dynamics are governed by a two- 
dimensional heat diffusion equation. Based on the concept of maximum likeli- 
hood estimate, the problem was converted into the optimal control problem. A 
difficulty arises in the computational burden involved in computing the gradi- 
ent of the associated likelihood ratio functional. In this paper, to avoid 
this difficulty, the steady state model was introduced. The validity of the 
estimator proposed here was demonstrated through numerical experiments. 
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Appendix 1: The Proof of Theorea 1 


hot ^ ^ be a sequence such that 


VV * J t (e t ). 


From (3.8), it yields that, for 0^, 0 €0, 


Noting that 


and 


we obtain 




+ 7 I H[u(0 ) + u(0 ) ] I . N ) 
1 n R (m) 


* UH[u(0 ) - u(0 ) ) I , 

n R (m) 


I u(0 )l , < C, 

n n l (G)~ 1 


E{-^- ly( t)l £ C£ for a fixed t € T, 


E{ KG ) ~ J,(0)|} < C, 1 u(0 ) - u(0)l , , 

j n h^(G) 


Using the compactness method, it can be shown that, under the condition 1 


u(0n) + u(0) strongly in H^(G). 


(A. 1) 

(A.2) 

(A. 3) 
(A. 4) 

(A. 5) 
C-l), 
(A. 6) 


Hence, we have 
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E{ |J t (0 n ) - J t (0)|} > 0. 


(A. 7 ) 


From (A. 1 ) and (A. 5), we may find 0 = 0 t * The proof has been completed. 
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Appendix 2. The Proof of Theorem 2 


A 

if r ho >pt. imjl solution 0 exists, the following variational inequal- 


ity nolds . 


7 9 J t (0 t ),(e " 9 t ) 2 0 


(A. 8) 


where V denotes the gradient of J (9 ) . In order to obtain V.J (9), 

^ t 0 t 

let us derive the local variation of cost functional (3.8). Let 0 A be the 
perturbed parameter of 0 such that 


9 (C) ■ 6(C) + A 69(5) for £ € 3G, 


(A. 9) 


where A is a positive constant and 69(£) is a given real valued func- 
tion. Set as 


«J t (0) = J t (0 A ) - J t (9), 


(A. 10) 


6u(x;9) = u(x;0 A ) - u(x;0). 


(A. 1 l ) 


Applying (3.8) to the mean value theorem, (A. 10) becomes 


SJ t (9) = -H [Su(9 ) ] ' - H[u(0 )]} 


(A. 12) 


On the other hand, <5u(x;0) Is a solution of 


5u(x) =0 in G 


(A# 13a) 
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with 


8(0 3 ~ ~ Tn ' ~ ' " + ° - 0 ^))6u(U) 


*68(5) 3u(5) _ 

1 -8(C) - *68(0 3n on d °U 


3(6u(Q) 

3n 


0 on 3G • 


Introduce the so-called "adjoint system" of which state p(x;0) is 
tion of 

Ap(x) * 0 in G 

with 

9(5 ) O " 9(5 ))p(5 ) = 0 on 

= H*f^- - H[u(0)]] on 3G 
o n t K 


where H [• J denotes the adjoint operator of H [ • ] • By virtue of 

formula, from (A# 13) and (A. 14), we have 


H[6u(8)r {—i " H fu(0)]} 


X [ c 59(0 {e(5)-l){9(5)+A«e(5)-l}- , 


x ’ 9 ) 3p(C ;8 ) 

3n 3n 


From (A. 10) and (A. 13), 5 3^(9) Is computed by 


(A. 13b) 
(A. 13c) 
a solu- 
(A. 14a) 
(A. 14b) 
(A. 14c) 
Green's 

(A. 15) 
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<5 J t (0 ) = 


^ L 69( ° {9(O~-lH0(e)+X60(O-l} 


. 3u(c;0) 3 p(C; 0) Jr 

X 9n 3n * 


(A. 16) 


Thus, we obtain 

V A J O).60 = 11m i 5J (0) 

y t , . a c 

x+o 


= - / 69(C) j 

a Gy {0(0-1} 

9u(C ;0) 9p(C;0) Jr 

~T^ ^ — dc - 


(A. 17) 


By setting 9 * 9 and 60 = 0 - 0^ in (A. 17), we can obtain the varia- 
tional inequality (3.11). 
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cost functional with respect to the perturbation of parameters, the estimation 
mechanism Ls proposed in a form of recursive computations. Finally, the feasi- 
bility of the estimator proposed here is demonstrated through results of digital 
simulation experiments. 
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